home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / graphic / tpega.zip / SURFACE.PAS < prev    next >
Pascal/Delphi Source File  |  1986-02-05  |  7KB  |  267 lines
  1. program surface;
  2. {taken from page 184 PC tech Journal Vol. 3 No. 11}
  3. var
  4.   xctr, yctr : real;
  5.  
  6. { *** NOTE: The following definitions are part of my math library.  They may or
  7.   may not be needed, depending on which function you put into function f(x,y).
  8.   You may want to omit these lines if you have your own functions library to
  9.   put in here.  For the sake of convenience, it would be wise to put these
  10.   definitions in a file and include it at compile time with the i directive. *}
  11.  
  12. { This program has been modified for the EGA by Kent Cedola, if you like a
  13.   copy of all of the EGA graphic routines for Turbo Pascal, send Disk and SASE
  14.   to Kent Cedola, 2015 Meadow Lake Court, Norfolk, VA 23518   * FREE *        }
  15.  
  16. const
  17.     Tiny: real = 1E-38;
  18.  
  19. {$I GPPARMS.P }
  20. {$I GPINIT.P  }
  21. {$I GPTERM.P  }
  22. {$I GPCOLOR.P }
  23. {$I GPMOVE.P  }
  24. {$I GPPLOT.P  }
  25. {$I GPLINE.P  }
  26.  
  27. function ArcSin(trig : real) : real;
  28. begin
  29.   ArcSin := ArcTan(trig / sqrt(1 - sqr(trig)));
  30. end; { function ArcSin(trig : real)  }
  31.  
  32. function ArcCos(trig : real) : real;
  33. begin
  34.   ArcCos := ArcTan(Sqrt(1 - Sqr(trig)) / trig);
  35. end; { function ArcCos(trig : real) }
  36.  
  37. function Tan(trig : real) : real;
  38. begin
  39.   Tan := Sin(trig) / Cos(trig);
  40. end;
  41.  
  42. function E(x,y : real) : real;   {raises x to the y pwr -- x^y in basic}
  43. begin
  44.   E := Exp(y * Ln(x));
  45. end; { function E(x,y }
  46.  
  47. { *** END OF MATH LIBRARY ***}
  48.  
  49. {Define function to be graphed - z=f(x,y)}
  50. function f(x,y:real): real;
  51. var
  52.   x1, y1 : real;
  53. begin
  54.   x1:=abs(x-xctr);  y1:=abs(y-yctr);
  55. { Function:  Enter your own below.  Notice that all others must be commented out
  56.   1. Pond Ripples
  57.   2. Stretched Planar Curve
  58.   3. Rippled Plane
  59.   4. Inverse Cone
  60.   5. Egg Carton
  61.   6. Pond Ripples 2
  62.   f := sin(sqrt(x1*x1+y1*y1));
  63.   f := sqrt(100-sqrt(x1*x1+y1*y1));
  64.   f := exp(-(x*y+y*y)/90)*cos((x*x+y*y)/40);
  65.   f := sqrt(x1*x1+y1*y1);                                   }
  66.   f := sin(x)+cos(y);                                       {
  67.   f := Sin(sqrt(x1*x1+y1*y1))/(sqrt(x1*x1+y1*y1)+tiny);
  68.   f := sin(1/x+1/(x*y));}
  69. end; { function f(x,y:real) }
  70.  
  71.  
  72. Const
  73.   xdiv = 40;            {number of subdivisions of each axis}
  74.   ydiv = 60;
  75.   xeye = 100;           {eye position}
  76.   yeye = 3;            {xeye and yeye should be positive}
  77.   zeye = 8;
  78.  
  79. var
  80.   i,j                  : integer;
  81.   xmin, xmax,
  82.   ymin, ymax,
  83.   zmin, zmax,
  84.   xdif, ydif, zdif     : real;
  85.   p, q                 : array[0..xdiv,0..ydiv] of integer;
  86.   y, z                 : array[0..xdiv,0..ydiv] of real;
  87.  
  88. {input extreme values for x and y}
  89.  
  90. procedure Input_Domain;
  91. begin
  92.   write('Enter smallest value of x  '); readln(xmin);
  93.   write('Enter largest value of x  '); readln(xmax); xdif:= xmax - xmin;
  94.   Write('Enter smallest value of y  '); readln(ymin);
  95.   write('Enter largest value of y  '); readln(ymax); ydif:= ymax - ymin;
  96.   xctr:= xmin + xdif/2;        yctr:= ymin + ydif/2;
  97. end; { procedure Input_Domain }
  98.  
  99. {evaluate function at grid points and project to view plane}
  100.  
  101. procedure Evaluate_and_Project;
  102. var
  103.   xtemp, xtemp1, xtemp2,
  104.   ytemp, ytemp1,
  105.   ztemp,
  106.   xavg,  yavg              : real;
  107.  
  108. begin
  109.  xavg:= (xmax+xmin)/2;  yavg:=(ymax+ymin)/2;
  110.  for i:=0 to xdiv do
  111.    for j:=0 to ydiv do begin
  112.      xtemp:=xmin+i*xdif/xdiv;
  113.      ytemp:=ymin+j*ydif/ydiv;
  114.      ztemp:=f(xtemp,ytemp);
  115.      xtemp1:=xeye - xtemp;
  116.      ytemp1:=yeye - ytemp;
  117.      y[i,j]:= (xeye - xavg)*(xeye*ytemp - yeye*xtemp) /
  118.              ((xeye - xavg)*xtemp1+(yeye-yavg)*ytemp1);
  119.      if y[i,j] <> yeye then
  120.        z[i,j] := zeye + (zeye - ztemp)*(Y[i,j] - yeye)/ytemp1
  121.      else begin
  122.        xtemp2 := yeye*(yavg - yeye) / (xeye - xavg);
  123.        z[i,j] := zeye + (zeye - ztemp)*(xtemp2 - xeye) / xtemp1;
  124.      end;
  125.    end;
  126. end; { procedure Evaluate_and_Project }
  127.  
  128. {determine projected extrema}
  129.  
  130. procedure Find_Extrema;
  131. var
  132.   ytemp, ztemp  : real;
  133. begin
  134.   ymax:= y[0,0];  ymin:=ymax;
  135.   zmax:= z[0,0]; zmin:=zmax;
  136.   for i:=0 to xdiv do
  137.     for j:=0 to ydiv do begin
  138.       ytemp:=y[i,j];
  139.       ztemp:=z[i,j];
  140.       if ytemp>ymax then ymax:=ytemp;
  141.       if ytemp<ymin then ymin:=ytemp;
  142.       if ztemp>zmax then zmax:=ztemp;
  143.       if ztemp<zmin then zmin:=ztemp;
  144.     end;
  145. end; { procedure Find_Extrema }
  146.  
  147. procedure Scale_to_Screen;
  148. var
  149.   dy, dz  : real;
  150. begin
  151.   dy:=(Ymax-ymin)/639; dz:=(Zmax - zmin)/349;
  152.   for i:=0 to xdiv do
  153.     for j:=0 to ydiv do begin
  154.       p[i,j]:=round((y[i,j] - ymin) / dy);
  155.       q[i,j]:=349 - round((z[i,j] - zmin) / dz);
  156.     end;
  157. end; { procedure Scale_to_Screen }
  158.  
  159. {exchange coordinates of two points}
  160.  
  161. procedure Swap(var x1,y1,x2,y2:integer);
  162. var
  163.   temp:integer;
  164. begin
  165.   temp:=x1; x1:=x2; x2:=temp;
  166.   temp:=y1; y1:=y2; y2:=temp;
  167. end; { procedure Swap(var x1,y1,x2,y2:integer) }
  168.  
  169. {draws blank horizontal line}
  170.  
  171. procedure Line(x0,x1,y:integer);
  172. begin
  173.   GPCOLOR(Black);
  174.   GPMOVE(x0,y);
  175.   GPLINE(x1,y);
  176. end; { procedure Line(x0,x1,y:integer) }
  177.  
  178. {blanks triangle}
  179.  
  180. procedure Triblank(x0,y0,x1,y1,x2,y2:integer);
  181. var
  182.   x3, x4,
  183.   dx1, dx2, dy1, dy2,
  184.   inc1, inc2,
  185.   nx1, nx2              : integer;
  186.  
  187.   procedure Blank(y:integer);
  188.   begin
  189.     while y0<y do begin
  190.       nx1:=nx1+dx1;
  191.       if nx1>dy1 then
  192.         repeat
  193.           x3:=x3+inc1;
  194.           nx1:=nx1-dy1;
  195.         until nx1<=dy1;
  196.       nx2:=nx2+dx2;
  197.       if nx2>dy2 then
  198.         repeat
  199.           x4:=x4+inc2;
  200.           nx2:=nx2 - dy2;
  201.         until nx2<=dy2;
  202.       y0:=y0+1;
  203.       line(x3,x4,y0);
  204.     end;
  205.   end; { procedure Blank(y:integer) }
  206.  
  207. begin
  208.   if y1<y0 then swap(x0,y0,x1,y1);
  209.   if y2<y0 then swap(x0,y0,x2,y2);
  210.   if y2<y1 then swap(x1,y1,x2,y2);
  211.   dy1:=y1-y0;    dy2:=y2-y0;
  212.   if x1<x0 then inc1:=-1 else inc1:=1;
  213.   if x2<x0 then inc2:=-1 else inc2:=1;
  214.   dx1:=abs(x1-x0);   dx2:=abs(x2-x0);
  215.   x3:=x0;  x4:=x0;
  216.   nx1:=dy1 div 2;    nx2:=dy2 div 2;
  217.   blank(y1);
  218.   if x2<x1 then inc1:=-1 else inc1:=1;
  219.   x3:=x1;  dy1:=y2-y1;
  220.   dx1:=abs(x1-x2);   nx1:=dy1 div 2;
  221.   blank(y2);
  222. end; { procedure Triblank(x0,y0,x1,y1,x2,y2:integer) }
  223.  
  224. {Draws box with blank interior}
  225.  
  226. procedure DrawBox(x1,y1,x2,y2,x3,y3,x4,y4 : integer);
  227. begin
  228.   triblank(x1,y1,x2,y2,x3,y3);
  229.   triblank(x2,y2,x3,y3,x4,y4);
  230.   GPCOLOR(Green);
  231.   GPMOVE(x1,y1); GPLINE(x2,y2);
  232.   GPMOVE(x1,y1); GPLINE(x3,y3);
  233.   GPMOVE(x2,y2); GPLINE(x4,y4);
  234.   GPMOVE(x3,y3); GPLINE(x4,y4);
  235. end; { procedure DrawBox(x1,y1,x2,y2,x3,y3,x4,y4 : integer) }
  236.  
  237. {Draws surface}
  238.  
  239. procedure Graph;
  240. var
  241.   x1,x2,x3,x4,y1,y2,y3,y4 : integer;
  242. begin
  243.   GPINIT;
  244.   for i:=0 to xdiv-1 do
  245.     for j:=0 to ydiv-1 do begin
  246.       x1:=p[i,j];    x2:=p[i+1,j];
  247.       x3:=p[i,j+1];  x4:=p[i+1,j+1];
  248.       y1:=q[i,j];    y2:=q[i+1,j];
  249.       y3:=q[i,j+1];  y4:=q[i+1,j+1];
  250.       drawBox(x1,y1,x2,y2,x3,y3,x4,y4);
  251.     end;
  252. end; { procedure Graph }
  253.  
  254.  
  255. begin
  256.   GPPARMS;
  257.   GPINIT;
  258.   GPCOLOR(Green);
  259.   input_Domain;
  260.   Evaluate_and_Project;
  261.   Find_Extrema;
  262.   Scale_to_Screen;
  263.   Graph;
  264.   repeat until keypressed;
  265.   GPTERM;
  266. end.
  267.